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Abstract. A uniform strategy to derive metric tensors in two spatial dimension 
for interpolation errors and their gradients in L p norm is presented. It generates 
anisotropic adaptive meshes as quasi-uniform ones in corresponding metric space, with 
the metric tensor being computed based on a posteriori error estimates in different 
norms. Numerical results show that the corresponding convergence rates are always 
optimal. 
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1 Introduction 

Generation of adaptive meshes is now the standard option in most software packages. 
Traditionally, isotropic mesh adaptation has received much attention, where regular mesh 
elements are only adjusted in size based on an error estimate. However, for problems 
with anisotropic solutions (with, say, sharp boundary or internal layers), the shape of 
elements can be further optimized and an equidistribution of a scalar error density is 
not sufficient to ensure that a mesh is optimally efficient [14] . Indeed anisotropic meshes 
have been used successfully in many areas, for example in singular perturbation and flow 
problems @H51 EH 1221133031 and in adaptive procedures [21 El HOI HU [23l [35l [36] . For 
anisotropic mesh adaptation, the common practice is to generate the needed anisotropic 
mesh as a quasi-uniform one in the metric space determined by a tensor (or a matrix- valued 
function), always called monitor function or metric tensor. Both the monitor function 
(denoted by the letter M) and metric tensor (denoted by the calligraphy letter M) play the 
same role in mesh generation, i.e., they are used to specify the size, shape, and orientation 
of mesh elements throughout the physical domain. The only difference lies in the way 
they specify the size of elements. Indeed, the former specifies the element size through the 
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equidistribution condition, while the latter determines the element size through the unitary 
volume requirement. Readers could regard the metric tensor as normalization for the 
monitor function. Examples of anisotropic meshing strategies include blue refinement \29\ 
[30] . directional refinement [36], Delaunay-type triangulation method [7,8, 1TTT35] . advancing 
front method [19], bubble packing method [H], local refinement and modification 
variational methods [21IT7JIM1E3I2H1EI] , and so on. Readers are referred to [T5] and 
for an overview. 

Among these meshing strategies, the definition of the metric tensor (or monitor function) 
based on the Hessian of the solution seems widespread in the meshing community [TjllH- 
[El[20l[2ll[23H26l[371[39]. Especially, Huang and Russell [26] propose the monitor function 



M = det ( I+-\H(u) 
a 



d+p(2-m) 
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I+-\H(u)\ d+pi2 - m) I+-\H(u)\ , (1.1) 
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for the interpolation error in W m ' p norm (m = 0, 1, p G [1, +oo)), where d stands for the 
spatial dimension. Set % = I + MH(u)\, when d = 2, 

1 mp 

M m>p = det(%)~ ^+p(2-m) ||^||2+P(2- m )^. (1.2) 

Separately, it becomes 

m , p = det(ny^^n, (1.3) 

for the interpolation error in L p norm and 

Mi iP = det(^)~^ \\U\\^U. (1.4) 

for the gradient of interpolation error in L p norm. 

The objective of this paper is to give a unified strategy deriving metric tensors in two 
spatial dimension for interpolation error and its gradients in L p norm. The development 
begin with the error estimates [32] for I? norm and our recent work [42J for H 1 norm on 
linear interpolation for quadratic functions on triangles. These estimates are anisotropic 
in the sense that they allow a full control of the shape of elements when used within a 
mesh generation strategy. Using the relationship between different norms, a posterior error 
estimates for other norms (W m ' p ,m = 0, 1, p / 2) can be gained. We will apply these 
error estimates to formulate corresponding metric tensors in a unified way. The procedure 
is based on two considerations: on the one hand the anisotropic mesh is generated as 
a quasi-uniform mesh in the metric tensor. On the other hand, the anisotropic mesh is 
required to minimize the error for a given number of triangles. To compare with those 
existing methods, we list our main results using monitor function style, that is 

M o yx) = det(H)~W^n, (1.5) 

for interpolation errors in L p norm and 

M£ p (x) = det(H)~^tT(H)^H, (1.6) 
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for gradient of interpolation errors in L p norm. To sum up, the metric tensor can be 
expressed by 

1 mp 

M" jp (x) = det(ft) - 2 +f( 2 - m > tr(ft) 2+ P (2- m j 7£ j (1.7) 

for the W m,p norm (m = 0, 1, p £ (0, +oo]) of the interpolation error. 

The paper is organized as follows. In Section 2, we describe the anisotropic error 
estimates on linear interpolation for quadratic functions on triangles obtained in our recent 
work |42j . The formulation of the monitor function and metric tensor is developed in 
Section 3. Numerical results are presented in Section 4 to illustrating our analysis. Finally, 
conclusions are drawn in Section 5. 



2 Estimates for interpolation error and its gradient 

As we know, the interpolation error depends on the solution, the size and shape of the 
elements in the mesh. Understanding this relation is crucial for the generating efficient 
meshes for the finite element method. In the mesh generation community, this relation 
is studied more closely for the model problem of interpolating quadratic functions. This 
treatment yields a reliable and efficient estimator of the interpolation error for general 
functions provided a saturation assumption is valid [51[TS]. For instance, Nadler [32] 
derived an exact expression for the L 2 -norm of the linear interpolation error in terms of 
the three sides £±, £ 2 , and £3 of the triangle K, 

\\ u ~ u i\\h(K) = (di + d 2 + d^j + did 2 + d 2 d 3 + did 3 , (2.1) 

where \K\ is the area of the triangle, di = £%• H£% with H being the Hessian of u. Assuming 
u = \±x 2 + \ 2 y 2 , D'Azevedo and Simpson |13| derived the exact formula for the maximum 
norm of the interpolation error 

n , s,|2 Di 2 D 23 D 31 

where Dj, = ^•diag(Ai, X 2 )£j- Based on the geometric interpretation of this formula, they 
proved that for a fixed area the optimal triangle, which produces the smallest maximum 
interpolation error, is the one obtained by compressing an equilateral triangle by factors 
\/Ai and y/X^ along the two eigenvectors of the Hessian of u. Furthermore, the optimal 
incidence for a given set of interpolation points is the Delaunay triangulation based on 
the stretching map (by factors vA~i an d VA2 along the two eigenvector directions) of the 
grid points. Rippa [38] showed that the mesh obtained in this way is also optimal for the 
L p -norm of the error for any 1 < p < 00. 

The element-wise error estimates in the following theorem are developed in |42| using 
the theory of interpolation and proper numerical quadrature formula. 
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Theorem 2.1. Let u be a quadratic function anduj is the Lagrangian linear finite element 
interpolation of u. The following relationship holds: 

1 3 

l|V(« - uMl^K) = 7^7[ • H£ i+ 2) 2 m 2 , (2.3) 



48 K 

i=i 

where we prescribe i + 3 = i,i — 3 = i. 

To get the a posteriori error estimate of the interpolation error in LP and W 1,p norms 
for p 7^ 2, we need some lemmas below. 

Lemma 2.1. For any d positive numbers a±, ■ ■ ■ , a^, the inequalities 

(£«?)* s-i-i (£«$)*■ (2 - 4) 

i=i i=i j=i 

and 

^ 4 (£<) 4 £(£#£(£<•?)' < 2 - 5) 

j=i j=i j=i 

ao/d /or numbers < p < 2 and p > 2, respectively. 

Proof. We just give the proof for the case < p < 2, it is similar for the case p > 2. 
For any number < p < 2, 

holds due to the Jensen's inequality. From the generalized arithmetic-mean geometric- 
mean inequality, for any positive numbers a\, ■ ■ ■ , aa, 

d ^ i d 



a i I a i I 



i=i i=i 

Then 



□ 



To sum up, for any d positive numbers a\, ■ ■ ■ , ad, the inequalities 

^(£s 2 ) f <(£«^<^(£«^ (^) 



j=i j=i j=i 



holds for any numbers p > 0, where C p = 1 for < p < 2 and dp 2 for p > 2, C p = d?> 
for < p < 2 and 1 for p > 2. 



i i 
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Lemma 2.2. J3]/ For any p £ (0, +00] and any non-negative v G P2(K) it holds 



with 



C,, p \K\v 1 \\v\ 
VP 



C p = l 



LHK) 



< V 



Lp(K) 



< CJK\ p \\v\\ti 



L 1 (K) 



if0<p<l, 



C p = (d + l)(d + 2)(d\)p(]J(P + 3)) " ifl<P<+oo, 
Coo= lim C p = (d+l){d + 2), 

p— >+oo 

C1/00 = hm Ci/ p = 1. 



(2.7) 



2.1 Estimates for interpolation errors in LP norm 

We consider the error of linear interpolation e = u — uj for a quadratic function u on K. 
Since the function e is quadratic on K, we can apply Lemma 12.21 to obtain 



a 



-Vp\ 
i/p I 



K\ 



< llel 



LP(K) 



<C p \K\p \\e\\ L i {K) . 



(2i 



Set p = 2, 



or 



l#1 2 IMlLi(iO < \\4lz(K) < C 2 \K\ 2 \\e\\ L i {K) , 



C 2 2 \\ e \\L*(K) < \\4l\k) < \K\ 2 IMIl 2 (^)- 



Combine (f2T%j) and ([23]), we get 

^i7 P /P C2 _1 |^l^^l|e||L2(K) < ||e||LP(jc) < C p \K\p-^\\e\\ L 2 {K) . 



(2.9) 



(2.10) 



In this article, A ~ i3 stands for that there exist two constants C_ and C such that 

CA < B < CA, 

where the two constants C_ and C may depend on the prescribed error, the index p, the 
dimension d, and the numbers of elements N, however are independent of function at 
hand. So (|2,10p can be rewritten as 



\ e \\LP(K) 



1_ 1 

\K\p 2 ||e| 



Together with the expression ()2.ip for the L 2 norm of the linear interpolation error derived 
by Nadler [32], we have the a posteriori error estimate in LP norms as follows: 



\Lp(K) 



\ K V 1 \\4h[ K ) 



2 

\K\p 
180 



6 

[(£■ 

i=i 



di) + d\d 2 + d 2 d 3 + d\d 3 



. (2.11) 
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2.2 Estimates for gradient of interpolation errors in LP norm 



Now we consider the gradient of linear interpolation error Ve = V(u — uj) for a quadratic 
function it. Since the function 

de \2 



is quadratic on K, we can apply Lemma 12.21 to obtain 



1 1/2 ■> 

1 1 LP/ 2 (K) ~ °2/p 



1 1/2 

u 3l\L^(K) ~ ^2/p 



L\K) 



(2.12) 



and 



1 1/2 



< cH*\K\p~*\\v 



1 1/2 



de 



dxj 



L 2 (K) 



(2.13) 



Since 



IVpII p 



i=i ' 





de 
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d Xj 



'*=£ 



i=i 



||P/2 

u illLP/2(i<r)' 



then together with (|2.12p and (j2. 13|) . we have 

d 



t-l/P| 
y 2/p I 



i=i 



1 ^ e 


P X 




L 2 {K)J 



< HVel 



^1/2, 



J'=l 



II ^ e 


P X 




L 2 (K)J 



From (|2.6p . the inequality 



-P^2/p 

holds, or simply 



C.C-^lXl^^Vell^^) < HVell^^) < C^J^K^-^e^^ 



|Ve| 



LP{K) 



1_ l 

~ 2 llVel 



Together with the a posteriori error estimate (|2.3p of the interpolation error in H 1 (= W 1,2 ) 
norm, we have the a posteriori error estimate in W 1 '? norms as follows: 



|Ve| 



LP(K) 



LK1* llVe 



\h(K) 



\K\^ 1 ^Y.(i l+ i-H K e l+2 ) 2 \e 

i=l 



2\n |2 



l^|--2 3 
A p 
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i+2, 



2\n |2 



(2.14) 
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3 Metric tensors for anisotropic mesh adaptation 



We now use the results of Section [2] to develop metric tensors for interpolation errors and 
their gradients in LP norm in a unified way. As a common practice in anisotropic mesh 
generation, the metric tensor, A4(x), is used in a meshing strategy in such a way that 
an anisotropic mesh is generated as a quasi-uniform mesh in the metric space determined 
by .M(x). Mathematically, this can be interpreted as the shape, size and equidistribution 
requirements as follows. 

The shape requirement. The elements of the new mesh, Th, are (or are close to 
being) equilateral in the metric. 

The size requirement. The elements of the new mesh Th have a unitary volume in 
the metric, i.e., 

/ ^det(M(x))dx =1, \/K G %. (3.1) 
Jk 

The equidistribution requirement. The anisotropic mesh is required to minimize 
the error for a given number of mesh points (or equidistribute the error on every element). 

Notice that to derive the monitor function, we just need the shape and equidistribution 
requirements. 

3.1 Metric tensors for gradients of interpolation errors in L p norm 




Figure 1: Affine map x = Tkx from K to the reference triangle K. 

We derive the monitor function M(x) first. Assume H(u) be a symmetric positive definite 
matrix on every point x, this assumption will be dropped later. Set M(x) = C(x)H(u). 
Consider the L? projection of H(u) on K, denoted by Hk, then so does Mk- Since Hk is 
a symmetric positive definite matrix, we consider the singular value decomposition Hk = 
R T AR, where A = diag(Ai,A2) is the diagonal matrix of the corresponding eigenvalues 
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(Ai, A2 > 0) and R is the orthogonal matrix having as rows the eigenvectors of Hk- Denote 
by Fk and tx the matrix and the vector defining the invertible affine map x = .F^ (x) = 

-Fr-x + tx from the generic element K to the reference triangle K (see Figure [1]). 

1 x 

Obviously, M K = C K H K . Let M K = F%F K , then F K = C^A^R. Mathematically, the 
shape requirement can be expressed as 



\£i\ = L and cos 9i 



L 2 ~ 2' ' ' ' 



(3.2) 



where L is a constant for every element K. Enforcing the shape requirement, we get 

2 o Q 

\K\p~ 



|Ve| 



Lp(K) 
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i=l 



-2 3 



-2 3 



Notice that, 



\K\ 



\K\ 



C Ky /det(H K y 



we have 



then 



|Ve| 



LP W '" 192Cf 
L 4 |i^|i~ 2 det(ff/<) 
192C K p 



i=i 



L 4 |if|p det^) p tr(gg) 

det(ilK) 



192C K p 



L 4 |iq? det(tf^) vtr{H K ) 



1+* 
192C^ p 



det(il x ) ftr(^) 



IVpII p 



|Ve| 



LP (A") 



)P/2 



det(g^)-2tr(g^) 
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To satisfy the equidistribution requirement, let 



Veil? 



LP(K) 



(£ek)/j\r = «w 



KeT h 



where N is the number of elements of Th- Then 

C K ~ det(H K y^tr(H K )^. 

So M(x) could be the form 

M(x) = det(H)~^tr(H)^+*H(u), 

since M(x) can be modified by multiplying a constant. Since it corresponds the gradient 
of interpolation errors in LP norm, we denote it by M" p (x). 

To establish the metric tensor .M™ p (x), set .M" p (x) = 0i jP M" (x), at this time, the 
size requirement (|3.ip should be used, which leads to 








where 



PliP (x) = ^det(M^(x)) 



Summing the above equation over all the elements of Th, one gets 



6i,p&i,p — N, 



where 




Thus, we get 
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N 



0\, P 



and as a consequence, 



IS 1 v 
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3.2 Metric tensor for the interpolation errors in L p norm 

Using the error estimates ()2.1ip for interpolation errors in L p norm and the shape require- 
ment (|3.2p . we have 



|e" 2 



, 2 3 



i=l 

3 



n i=i i=i 
L 4 \K\p L*\K\p 1 



1 ^n 2 2+- i 2+- 

l5C K p det{H)p C K p det{H)i 



Then, 



Wils i Nl"!/v)''' " ~ . .,,-1 . ' • 



|ef r „_ = (|| e ||? B ^) p/2 



_• 



C£ +1 det(tf) 
To satisfy the equidistribution requirement, let 

II< p( k) = (E4)av = ^. 

Using similar argument in last subsection, we easily get monitor functions 

M^ p ( X ) = det(Hy^H(u), 

and metric tensors 

M£ p (x) = —det(H)-W)H(u), 
a o, P 

for the interpolation errors in LP norm. 



3.3 Practice use of metric tensor 

So far we assume that H (u) is a symmetric positive definite matrix at every point. However 
this assumption doesn't hold in many cases. In order to obtain a symmetric positive 
definite matrix, the following procedure are often implemented. First, the Hessian H is 
modified into \H\ = R T diag(|Ai|, IA2D-R by taking the absolute value of its eigenvalues 
( [22] )• Since \ H\ is only semi-positive definite, Ai^ p cannot be directly applied to generate 
the anisotropic meshes. To avoid this difficulty, we regularize the expression with the 
flooring parameter a„ hP > (see, e.g., |24|). Replacing \H\ with 

W = a m , P I + \H\, 
we get the modified metric tensors, also denoted by .M" , that is 

]\[ 1 mp 

AC, P ( X ) = det(H) a+PCa-m)tr(W) a -t*P-m)ft f (3.3) 

which are suitable for practical mesh generation. 



10 



10* 
nbt 



1( 



10 4 
nbt 



io- 



ta) (b) 

Figure 2: Example 1: Interpolation error and its gradient in L p norm 

3.4 Comparison with existing methods using monitor function style 

When m = 0, the new monitor function (jl.5p is in fact the same with (jl.3p in |25U26j . 
Chen, Sun and Xu [12] proved that under suitable conditions, the error estimate 

ll n ~~ u i\\lp(q.) — CN~d II \J det H\\ V d , l<p<oo, 

1 

holds on the quasi-uniform mesh determined by the metric (det-ff) 2 p+<iH, where H is a 
majorant of the Hessian matrix, N is the number of elements in the triangulation and the 
constant C does not depend on u and N. This estimate is optimal in the sense that it is a 
lower bound if u is strictly convex or concave. Note that % can be chosen as a majorant 
of the Hessian matrix. 

When m = 1, the new monitor function Mf p (jl.6p is different with (|1.4p |26j that 
the former refers to tr(%) and the latter involves \\H\\. In some cases, the two monitor 
functions are pretty much alike. However, in other cases, the effect of the former is superior 
to the latter for mesh generation. Numerical results in [3D] have shown our approach's 
superiority for the error in H l norm. 

4 Numerical experiments 

In this section, we present some numerical results for three problems with given analytical 
solutions. The numerical results are performed by using the BAMG software [23]. Given 
a background mesh and an approximation solution, BAMG generates the mesh according 
to the metric tensor. The code allows the user to supply his/her own metric tensor defined 
on a background mesh. In our computation, the background mesh has been taken as the 
most recent mesh available. 
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(a) (b) 
Figure 3: Example 1: plots of the solution (a) and corresponding mesh (b) using -M^i 

Denote by nbt the number of triangles in the current mesh. The number of triangles is 
adjusted when necessary by trial and errors through the modification of the multiplicative 
coefficient of the metric tensors. 

Example 1 This example is to generate adaptive meshes for 

U(X) = ^TTCT n^ ' x E (0.1, 1) X (0.1, 1). (4.1) 

I _|_ e -200(- v /a:f+a;|-0.8) 

This function is anisotropic along the quarter circle x\ + x\ = 0.8 2 and changes sharply 
in the direction normal to this curve. A similar example was presented in [33] where the 
region is (0, 1) x (0, 1). In the current computation, each run is stopped after 15 iterations 
to guarantee that the adaptive procedure tends towards stability. We show in Figure [2] the 
LP norms of the interpolation error and its gradient using corresponding metric tensors, for 
p = 1, 2, 4, oo. For example, the curve p = 2 in (a) stands for the interpolation error using 
the metric tensor while p = oo in (b) stands for the gradient interpolation error 

using the metric tensor M\ )0 o- We see that the convergence rates for the interpolation 
error and its gradient are always nearly optimal, i.e. ||e||xp ~ N^ 1 and ||Ve||ij> ~ iV~ ' 5 . 
We also show in Figure [3] plots of the solution and corresponding mesh using the metric 
tensor Aii t i. 

Example 2 This example is to generate adaptive meshes for 

u(x) = x\x 2 +xl + tanh(10(sin(5a;2) - 2x 1 )), xG (-1,1) X (-1,1). (4.2) 

This function is anisotropic along the zigzag curve sin(5x2) — 2x\ = and changes sharply 
in the direction normal to this curve (taken from [3]). In the current computation, each 
run is stopped after 20 iterations to guarantee that the adaptive procedure tends towards 
stability. We show in Figure H] the LP norms of the interpolation error and its gradient 
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(a) (b) 
Figure 4: Example 2: Interpolation error and its gradient in LP norm 

using corresponding metric tensors, for p = 1,2,4, oo. As in Example 1, the convergence 
rates for the interpolation error and its gradient here are always nearly optimal. In Figure 
[5] we select 6 meshes with 4000 triangles generated by corresponding metric tensors. We 
can learn that the optimal meshes in different norms are different. For example, the mesh 
generated by the metric tensor Aii ;00 concentrates more triangles and nodes along the 
zigzag line. 

Example 3 (Taken from [30]) This example is to solve the boundary value problem of 
Poisson's equation 

-Au = /, x € ft = (-1.2, 1.2) x (-1.2,1.2), (4.3) 

with the Dirichlet boundary condition and the right-hand side term being chosen such 
that the exact solution is given by 

5 

E r x-\-y — ca x — y — dj n 

[(1 + e-E" 1 )- 1 + (1 + e-^Y 1 ] , (4.4) 

i=\ 

where a = 0,-0.6,0.6,-1.2,1.2; <& = 0,-0.6,0.6,-1.2,1.2. The solution exhibits ten 
sharp layers on lines x+y — Ci = and x — y — di = 0, i = 1, 2, • • • ,5, when e is small. In our 
computations, e is taken as 0.01. Numerical results in [40] have shown that our approach's 
superiority for the error in H 1 norm. In the current computation, each run is stopped after 
20 iterations to guarantee that the adaptive procedure tends towards stability, except that 
governed by A^i )00 , which need 30 iterations. We show in Figure [6] the LP norms of the 
interpolation error and its gradient using corresponding metric tensors, for p = 1,2,4, oo. 
As in Example 1 and Example 2, the convergence rates for the interpolation error and 
its gradient here are always nearly optimal. Another purpose to select this example is to 
describe the difference of finding layers using different norms. In Figure[7]we list meshes in 
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(d) (e) (f) 

Figure 5: Example 2: Meshes generated by the metric tensor M. m ,p for (a)m = 0,p = 1, 
(b)m = 0,p = 2, (c)m = 0,p = oo, (d)m = l,p = 1, (e)m = l,p = 2, (f)m = l,p = oo. 

different stage during one selected run governed by corresponding metric tensors. While 
in Figure [8] convergence history is shown. From the three figures we can learn that most of 
the metric tensors can quickly find the layers except the metric tensor Aix jOC when dealing 
with the complex problems, e.g., with multiple layers. 

5 Conclusions 

In the previous sections we have developed a uniform strategy to derive metric tensors in 
two spatial dimension for interpolation errors and their gradients in LP norm. The metric 
tensor for the LP norm of the interpolation error is similar to some existing methods. 
However, the metric tensor M™ p is essentially different with those metric tensors existed. 
There is a fine distinction between the new metric tensor Mi p and Mi p proposed by 
Huang and Russell [26] that the former refers to tr('H) and the latter involves \\7i\\. In 
some cases, the two metric tensors are pretty much alike. However, when dealing with 
the complex problems, e.g., with multiple layers, the effect of the former is superior to 
the latter for mesh generation. Numerical results show that the corresponding convergent 
rates are always almost optimal. 
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(a) (b) 
Figure 6: Example 3: Interpolation error and its gradient in L p norm 
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